Predicting range shifts of the giant pandas under future climate and land use scenarios

Abstract Understanding and predicting how species will respond to global environmental change (i.e., climate and land use change) is essential to efficiently inform conservation and management strategies for authorities and managers. Here, we assessed the combined effect of future climate and land use change on the potential range shifts of the giant pandas (Ailuropoda melanoleuca) in Sichuan Province, China. We used species distribution models (SDMs) to forecast range shifts of the giant pandas by the 2050s and 2070s under four combined climate and land use change scenarios. We also compared the differences in distributional changes of giant pandas among the five mountains in the study area. Our SDMs exhibited good model performance and were not overfitted, with a mean Boyce index of 0.960 ± 0.015 and a mean omission rate of 0.002 ± 0.003, and suggested that precipitation seasonality, annual mean temperature, the proportion of forest cover, and total annual precipitation are the most important factors in shaping the current distribution pattern of the giant pandas. Our projections of future species distribution also suggested a range expansion under an optimistic greenhouse gas emission, while suggesting a range contraction under a pessimistic greenhouse gas emission. Moreover, we found that there is considerable variation in the projected range change patterns among the five mountains in the study area. Especially, the suitable habitat of the giant panda is predicted to increase under all scenarios in the Minshan mountains, while is predicted to decrease under all scenarios in Daxiangling and Liangshan mountains, indicating the vulnerability of the giant pandas at low latitudes. Our findings highlight the importance of an integrated approach that combines climate and land use change to predict the future species distribution and the need for a spatial explicit consideration of the projected range change patterns of target species for guiding conservation and management strategies.


| INTRODUC TI ON
The earth is undergoing serious biodiversity loss due to humaninduced global environment change (Pimm et al., 2014;Tittensor et al., 2014), of which climate change and land use change rank among the most important direct drivers of such trends (Di Febbraro et al., 2019;Sala et al., 2000). This threat is predicted to become more intense in the near future due to accelerating global warming and intensifying habitat fragmentation (Sirami et al., 2017;Titeux et al., 2017). To avoid the further loss of biodiversity, effective conservation programs are required to mitigate the negative impact of these factors (Barnosky et al., 2011). A foundational role of biodiversity conservation research is to understand how these factors might affect the future distribution of species so as to efficiently inform conservation and management strategies for authorities and managers (Coreau et al., 2009;Maggini et al., 2014).
A common practice for assessing the impacts of climate and land use change on species is to use species distribution models (SDMs) to project the species range shifts under future environmental conditions (Long et al., 2021;Marshall et al., 2018;Prestele et al., 2021;Schweiger et al., 2012). Most studies have projected the future species distribution under the combination of "dynamic" (i.e., change with time periods of projection) climate variables and "static" (i.e., remain unchanged) land use variables. However, due to the ongoing and continuing land use change, the static land use variables may not fully represent future species habitat suitability (Martin et al., 2013), which could potentially lead to unrealistic projections of future species distribution (Marshall et al., 2018;Sirami et al., 2017;Titeux et al., 2017). Moreover, the use of static land use has been criticized for neglecting the effect of land use change on the future distribution of a species (Barbet-Massin et al., 2012). Therefore, an integrated approach that combing dynamic climate variables and dynamics land use variables is essential to project species' future distribution and assess the effect of land use change on species' distribution shifts. With the increase in available land use datasets at fine spatial scales under multiple future scenarios (Hurtt et al., 2011;Li et al., 2016), this type of integrated approach has been widely used to estimate species' future distribution changes for many taxonomic groups, such as plants (García-Valdés et al., 2015;Zhang et al., 2017), insects (Marshall et al., 2018;Prestele et al., 2021), and mammals (Ma et al., 2021;Zamora-Gutierrez et al., 2018).
As an iconic species and global symbol of conservation, the giant panda (Ailuropoda melanoleuca) has undergone pronounced humandriven range contractions over the past 3000 years (Zhao et al., 2013), and now only lives in six isolated mountain ranges in Sichuan, Shaanxi, and Gansu provinces in southcentral China: the Qinling mountains of southeastern Shanxi province, the Minshan mountains of the southern Gansu province, the northwestern part of Sichuan province, and the Qionglai, Xiaoiangling, Daxiangling, and Liangshan mountains of southwestern Shanxi province. Over the past decades, China has established 67 panda nature reserves to protect this species from intensifying human activities and environmental changes (State Forestry Administration, 2021). Besides, two national-level environmental protection projects: the Grain-to-Green Program (GTGP) and the Natural Forest Conversion Program (NFCP) were implemented to converse agricultural land with forested land, which also contributes to the recovery and improvement of panda habitats . Due to these conservation efforts, the giant pandas have recently been downlisted from Endangered to Vulnerable by the International Union for Conservation of Nature Red List (Swaisgood et al., 2016).
However, SDMs predicted that it will face serious risks from future climate change as a result of the dramatic loss of suitable habitats and intensifying habitat fragment (Li et al., 2015;Li et al., 2017;Shen et al., 2015;Songer et al., 2012;Wang et al., 2018). Moreover, previous studies have reported that giant pandas mainly inhabit primary forests (Hong et al., 2015;Zhang et al., 2011), and the knowledge regarding the effect of human-induced land use patterns on the distribution and habitat selection of the giant panda at fine scales is well established (Bai et al., 2020;Wei et al., 2018;Yang et al., 2020).
However, the effect of climate change and land use change, as well as their interactions, on a large-scale distribution of the giant panda has rarely been explored. Tang et al. (2020) were the first to include dynamic land-use changes on a large scale to assess the relative importance of climate change and land use change in determining the historical distribution patterns of the giant panda, finding that land use change could offset some of the negative effects of climate change on the giant pandas. Yet, to date, no studies have integrated climate change and land use change to project the future distribution of the giant panda at large scales.
To fill this gap, in this study we explored the combined effects of climate change and land use change on the future distribution changes in giant pandas. Our specific objectives are as follows: (1) identify the environmental requirements (i.e., ecological niche) of the giant pandas; (2) assess its habitat suitability under current Administration, 2021). In total, 3428 occurrence records for the giant pandas in our study area were obtained from this survey (Tang et al., 2020). To be consistent with the spatial resolution of our climate and land use variables (see below) and minimize the sampling bias effect in the occurrence records dataset, we created the same 1 × 1 km grid cells across the study area as the climatic and land use layers and all the occurrence records were overlaid onto these grid cells. After removing duplicate records within each gird cell, we obtained 2068 occurrence records to model ecological niches for the giant pandas.

| Species distribution modeling
The maximum entropy algorithm (MaxEnt; Phillips et al., 2006) was used to make current and future projections of potentially suitable habitats for giant pandas. We chose MaxEnt because of its superior performance to model species distribution using presenceonly data compared to other algorithms (Elith et al., 2006).
Moreover, recent studies suggested that the tuned MaxEnt models can perform comparably to ensemble SDMs (Hao et al., 2020;Low et al., 2021). To improve the performance of MaxEnt and avoid overfitting, following Jarvie et al. (2021), we ran 24 Maxent models based on all possible combinations of eight regularization multipliers (i.e., 0.5, 1, 1.5, 2, 2.5, 3, 3.5, and 4) and three feature class options (i.e., linear, linear/quadratic, and linear/quadratic/ product). The performance of each model was evaluated using a non-spatial fivefold cross-validation. We used the corrected Akaike information criterion (AICc) to select the best-performing model as this metric reflects both model goodness of fit and complexity (Muscarella et al., 2014). All the models were developed using the ENMeval package (Version 2.0.3; Kass et al., 2021) with the ENMevaluate function in the R platform (v. 4.1.3; http://cran.rproje ct.org).
The performance of the optimal model was evaluated using the Boyce index (Boyce et al., 2002). We chose the Boyce index due to its superior performance compared to the two commonly used metrics, i.e., the area under the receiver operating characteristic curve (AUC) and the true skill statistic (TSS), for both the metrics may present a problem when presence-only data are used (Jiménez-Valverde, 2012;Leroy et al., 2018;Lobo et al., 2008).
The Boyce index ranges from −1 to 1, with positive values indicating better model performance and negative values indicating performance no better than a random model (Hirzel et al., 2006).
In addition, to assess whether the optimal model is overfitting, we also calculated the test point omission rate based on the minimum training presence value (OR MTP ). This metric was threshold dependent and range from 0 (models that are not overfitted) to 1.0 (models that are overfitted; Peterson et al., 2011). The relative importance of the eight explanatory variables from the optimal model was determined by using the "permutation contribution,"  (Phillips et al., 2006). One habitat suitability was then produced using the optimal model for the current period and for each of the four different future scenarios (i.e., two RCP scenarios [RCP2.6 and RCP8.5] at two time periods [the 2050s and 2070s]), respectively. Finally, all these maps were converted into binary presence-absence maps by using the threshold the maximum sum of sensitivity and specificity, which has been frequently recommended .

| Model performance and variable contributions
The settings (regularization multiplier = 0.5 and feature combination = linear/quadratic/product) yielded the optimal model (ΔAICc = 0) ( Figure 1). The optimal models had an excellent predictive performance, with a mean Boyce index of 0.960 ± 0.015.
Besides, the mean OR MTP value is 0.002 ± 0.003, indicating that our models were not overfitted. Among the eight selected predictor variables, precipitation seasonality (BIO15) had the highest contribution to our model, followed by annual mean temperature (BIO1), the proportion of forest area (FL), and total annual precipitation (BIO12), while the remaining four variables contributed little to the distribution of the giant pandas ( Figure 2).

| Habitat suitability under current and future environmental conditions
The predicted potential suitable habitat area for the giant pandas in the study area under current and future environmental conditions is presented in  Figures 3 and 4). Specifically, under RCP 2.6, we predicted a net gain of 4.20% (1086 km 2 ) in currently suitable habitats based on the projected loss of 18.17% (4690 km 2 ) and a gain of

F I G U R E 1
The ΔAICc (the corrected Akaike information criterion) values for the MaxEnt models under a range of feature combinations and regularization multipliers. The settings (regularization multiplier = 0.5 and feature combination = LQP) yielded the best-performing model (ΔAICc = 0). L, linear feature; Q, quadratic feature; P, product feature.

F I G U R E 2
The mean variable importance of the eight selected environmental variables is included in the optimal MaxEnt model. BIO1: Annual mean temperature, BIO4: Temperature seasonality, BIO12: Total annual precipitation, BIO15: Precipitation seasonality, CL, FL, SL, and UGSL are the proportion of area covered by cropland, forest, shrubland, and urban green spaces in grid cells, respectively.

| DISCUSS ION
In this study, we provide the first empirical test of the combined impacts of climate and land use change on the future distribution changes of giant pandas. Projections suggest that the future changes in suitable habitats of the giant pandas largely depend on the greenhouse gas emission scenario, that is, the giant pandas will experience range expansion under RCP 2.6, but will experience range contraction under RCP 8.5 for the whole study area. Moreover, we found that there is considerable variation in the projected range change It has been reported that changes in land-use patterns have an important effect on the distribution and habitat selection of the giant pandas as these changes can have either a positive or negative effect on the habitat of the giant pandas, especially for the forest cover, a key natural resource for the giant pandas (Bai et al., 2020;Tang et al., 2020;Wei et al., 2018). Consistent with these studies, we found that land use variables are important drivers in determining the current distribution of the giant pandas in our study area.
More importantly, among these land use variables, the proportion of the forest cover contributes the most to our optimal MaxEnt models. As with other species (e.g., Marshall et al., 2018), these findings highlight the importance of incorporating land use variables into the SDMs with only climate variables.
Previous studies have suggested that the giant pandas have shifted and will continue to shift toward high altitude and/or latitude in response to the ongoing climate change, which would lead to a dramatic loss of suitable habitat for the giant pandas ranging from 30% to 85% in the future (Li et al., 2015;Li et al., 2017;Shen et al., 2015;Songer et al., 2012;Wang et al., 2018). Consistent with these previous studies, our integrated models that combine climate and land use change also predicted that about 72% of the current suitable habitat of the giant pandas was distributed at high latitudes (i.e., Minshan and Qionglai mountains), and this proportion will further increase in the future. Nonetheless, compared to these previous studies, our integrated models predicted less range contraction, even range expansion under the low greenhouse gas emission scenario. For example, Li et al. (2015) predicted that at least 52.9% of the current suitable habitat of the giant pandas will be lost under future climate change, while our integrated models predicted that the suitable habitat area will decrease by at most 25% in the future.
These findings suggest that conservation intervention to manage habitat could offset the negative effects of climate change on the future distribution of giant pandas to some extent. together with the narrow range of the giant pandas, could lead to a heavy loss of suitable habitat in these two mountains (Williams et al., 2007), which also highlights that the prior implementation of effective conservation programs in these two mountains is of great urgency and significance. Despite that, the Liangshan mountains have not been included in the recently established Giant Panda National Park (Xu et al., 2017).
Finally, a better understanding of how climatic and land-use changes will interact to influence future species distributions provides a platform for more informed conservation management strategies. Compared to climatic variables, land use variables such as forest cover are more under the control of managers and decision-makers at regional and local scales. Over the past decades, many conservation efforts, such as the establishment of protected areas and the implementation of NFCP and GTGP, have greatly contributed to the improvement of forest cover and accordingly contributed to the nascent recovery of the panda (Swaisgood et al., 2016. Similarly, our findings indicated that these conservation efforts would also have positive impacts on the future distribution of the giant pandas a the land use patterns would offset some of the negative effects of changing climate. Our findings also highlighted the need for a spatial explicit consideration of the projected range change patterns of target species when assessing the effect of climate and land use change on the future distribution of species as managers and decision-makers should be informed when and where the target species is likely at risk. However, due to methodological limitations such as study area, data, and algorithms, our analyses might lead to uncertainties regarding the projected species distribution (Préau et al., 2022). Therefore, caution is warranted in generalizing our findings to local regions to avoid the risk of mismanagement and implement conservation strategies to mitigate climate change by adapting land use management (Oliver & Morecroft, 2014). writing -review and editing (equal). Zenjun Zhang: Supervision (equal); writing -review and editing (equal).

ACK N OWLED G M ENTS
This study was supported by the National Natural Science University.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.